# -*- coding: utf-8 -*-
"""
Created on Mon Dec 20 21:23:38 2021

@author: LSY
"""

from __future__ import division
import numpy as np


Nx = 78
def ode_wwtp_scaled(x_scale,u,w=np.zeros(78)): #w=np.zeros(36)
    
    xs = np.loadtxt('x0_78.txt')  
    min_x = np.zeros(Nx) 
    max_x = np.zeros(Nx)  
    lb = 0.5
    ub = 1.5
    for i in range(Nx):
        min_x[i] = lb*xs[i]
        max_x[i] = ub*xs[i]   
    delta_x = max_x - min_x
    x = x_scale*delta_x + min_x
#    
   # Input  16 variables
    KLa5 = u[0]
    Qa = u[1]
    Q0 = u[2]
    Z0 = u[3:16]
    
    x1 = x[0] 
    x2 = x[1] 
    x3 = x[2] 
    x4 = x[3] 
    x5 = x[4] 
    x6 = x[5]
    x7 = x[6] 
    x8 = x[7] 
    x9 = x[8] 
    x10 = x[9] 
    x11 = x[10] 
    x12 = x[11]
    x13 = x[12] 
    x14 = x[13] 
    x15 = x[14] 
    x16 = x[15] 
    x17 = x[16]
    x18 = x[17] 
    x19 = x[18] 
    x20 = x[19] 
    x21 = x[20] 
    x22 = x[21]
    x23 = x[22] 
    x24 = x[23] 
    x25 = x[24] 
    x26 = x[25] 
    x27 = x[26] 
    x28 = x[27] 
    x29 = x[28] 
    x30 = x[29] 
    x31 = x[30] 
    x32 = x[31]
    x33 = x[32] 
    x34 = x[33] 
    x35 = x[34] 
    x36 = x[35] 
    x37 = x[36]
    x38 = x[37] 
    x39 = x[38] 
    x40 = x[39] 
    x41 = x[40] 
    x42 = x[41]
    x43 = x[42] 
    x44 = x[43] 
    x45 = x[44] 
    x46 = x[45] 
    x47 = x[46]
    x48 = x[47] 
    x49 = x[48] 
    x50 = x[49] 
    x51 = x[50] 
    x52 = x[51]
    x53 = x[52] 
    x54 = x[53] 
    x55 = x[54] 
    x56 = x[55] 
    x57 = x[56]
    x58 = x[57] 
    x59 = x[58] 
    x60 = x[59] 
    x61 = x[60] 
    x62 = x[61] 
    x63 = x[62] 
    x64 = x[63] 
    x65 = x[64] 
    x66 = x[65] 
    x67 = x[66] 
    x68 = x[67] 
    x69 = x[68] 
    x70 = x[69] 
    x71 = x[70] 
    x72 = x[71]
    x73 = x[72] 
    x74 = x[73] 
    x75 = x[74] 
    x76 = x[75] 
    x77 = x[76] 
    x78 = x[77]
    
    Z01 = Z0[0]
    Z02 = Z0[1]
    Z03 = Z0[2]
    Z04 = Z0[3]
    Z05 = Z0[4]
    Z06 = Z0[5]
    Z07 = Z0[6]
    Z08 = Z0[7]
    Z09 = Z0[8]
    Z010 = Z0[9]
    Z011 = Z0[10]
    Z012 = Z0[11]
    Z013 = Z0[12]
 ### SI, SS, XI, XS, XBH, XBA, XP, SO, SNO, SNH, SND, XND, SALK   
    dx_scaledt =  1./delta_x*[
            1/96*( (Q0*Z01)/1000 + (Qa*x53)/1000 - (x1*(Q0 + Qa + 18446))/1000),
            1/96*((Q0*Z02)/1000 + (Qa*x54)/1000 - (x2*(Q0 + Qa + 18446))/1000 + (3*x4*(x8/(x8 + 1/5) + (4*x9)/(25*(x9 + 1/2)*(x8 + 1/5))))/(x4/x5 + 1/10) - (400*x2*x5*x8)/(67*(x2 + 10)*(x8 + 1/5)) - (64*x2*x5*x9)/(67*(x2 + 10)*(x9 + 1/2)*(x8 + 1/5))),
            1/96*((Q0*Z03)/1000 + (Qa*x55)/1000 - (x3*(Q0 + Qa + 18446))/1000 + (9223*x55*x66)/(500*((3*x55)/4 + (3*x56)/4 + (3*x57)/4 + (3*x58)/4 + (3*x59)/4))),
            1/96*((69*x5)/250 + (23*x6)/500 + (Q0*Z04)/1000 + (Qa*x56)/1000 - (x4*(Q0 + Qa + 18446))/1000 + (9223*x56*x66)/(500*((3*x55)/4 + (3*x56)/4 + (3*x57)/4 + (3*x58)/4 + (3*x59)/4)) - (3*x4*(x8/(x8 + 1/5) + (4*x9)/(25*(x9 + 1/2)*(x8 + 1/5))))/(x4/x5 + 1/10)),
            1/96*((Q0*Z05)/1000 - (3*x5)/10 + (Qa*x57)/1000 - (x5*(Q0 + Qa + 18446))/1000 + (9223*x57*x66)/(500*((3*x55)/4 + (3*x56)/4 + (3*x57)/4 + (3*x58)/4 + (3*x59)/4)) + (4*x2*x5*x8)/((x2 + 10)*(x8 + 1/5)) + (16*x2*x5*x9)/(25*(x2 + 10)*(x9 + 1/2)*(x8 + 1/5))),
            1/96*((Q0*Z06)/1000 - x6/20 + (Qa*x58)/1000 - (x6*(Q0 + Qa + 18446))/1000 + (9223*x58*x66)/(500*((3*x55)/4 + (3*x56)/4 + (3*x57)/4 + (3*x58)/4 + (3*x59)/4)) + (x6*x8*x10)/(2*(x10 + 1)*(x8 + 2/5))),
            1/96*((3*x5)/125 + x6/250 + (Q0*Z07)/1000 + (Qa*x59)/1000 - (x7*(Q0 + Qa + 18446))/1000 + (9223*x59*x66)/(500*((3*x55)/4 + (3*x56)/4 + (3*x57)/4 + (3*x58)/4 + (3*x59)/4))),
            1/96*( (Q0*Z08)/1000 + (Qa*x60)/1000 - (x8*(Q0 + Qa + 18446))/1000 - (132*x2*x5*x8)/(67*(x2 + 10)*(x8 + 1/5)) - (433*x6*x8*x10)/(48*(x10 + 1)*(x8 + 2/5))),
            1/96*((Q0*Z09)/1000 + (Qa*x61)/1000 - (x9*(Q0 + Qa + 18446))/1000 + (25*x6*x8*x10)/(12*(x10 + 1)*(x8 + 2/5)) - (96*x2*x5*x9)/(871*(x2 + 10)*(x9 + 1/2)*(x8 + 1/5))),
            1/96*( (Q0*Z010)/1000 + (Qa*x62)/1000 + (x5*x11)/20 - (x10*(Q0 + Qa + 18446))/1000 - (8*x2*x5*x8)/(25*(x2 + 10)*(x8 + 1/5)) - (637*x6*x8*x10)/(300*(x10 + 1)*(x8 + 2/5)) - (32*x2*x5*x9)/(625*(x2 + 10)*(x9 + 1/2)*(x8 + 1/5))),
            1/96*( (Q0*Z011)/1000 + (Qa*x63)/1000 - (x5*x11)/20 - (x11*(Q0 + Qa + 18446))/1000 + (3*x12*(x8/(x8 + 1/5) + (4*x9)/(25*(x9 + 1/2)*(x8 + 1/5))))/(x4/x5 + 1/10)),
            1/96*((141*x5)/6250 + (47*x6)/12500 + (Q0*Z012)/1000 + (Qa*x64)/1000 - (x12*(Q0 + Qa + 18446))/1000 + (9223*x64*x66)/(500*((3*x55)/4 + (3*x56)/4 + (3*x57)/4 + (3*x58)/4 + (3*x59)/4)) - (3*x12*(x8/(x8 + 1/5) + (4*x9)/(25*(x9 + 1/2)*(x8 + 1/5))))/(x4/x5 + 1/10)),
            1/96*( (Q0*Z013)/1000 + (Qa*x65)/1000 + (x5*x11)/280 - (x13*(Q0 + Qa + 18446))/1000 - (4*x2*x5*x8)/(175*(x2 + 10)*(x8 + 1/5)) - (631*x6*x8*x10)/(2100*(x10 + 1)*(x8 + 2/5)) + (1518823277841921*x2*x5*x9)/(360287970189639680*(x2 + 10)*(x9 + 1/2)*(x8 + 1/5))),

            1/96*((x1*(Q0 + Qa + 18446))/1000 - (x14*(Q0 + Qa + 18446))/1000),
            1/96*((x2*(Q0 + Qa + 18446))/1000 - (x15*(Q0 + Qa + 18446))/1000 + (3*x17*(x21/(x21 + 1/5) + (4*x22)/(25*(x22 + 1/2)*(x21 + 1/5))))/(x17/x18 + 1/10) - (400*x15*x18*x21)/(67*(x15 + 10)*(x21 + 1/5)) - (64*x15*x18*x22)/(67*(x15 + 10)*(x22 + 1/2)*(x21 + 1/5))),
            1/96*((x3*(Q0 + Qa + 18446))/1000 - (x16*(Q0 + Qa + 18446))/1000),
            1/96*((69*x18)/250 + (23*x19)/500 + (x4*(Q0 + Qa + 18446))/1000 - (x17*(Q0 + Qa + 18446))/1000 - (3*x17*(x21/(x21 + 1/5) + (4*x22)/(25*(x22 + 1/2)*(x21 + 1/5))))/(x17/x18 + 1/10)),
            1/96*((x5*(Q0 + Qa + 18446))/1000 - (3*x18)/10 - (x18*(Q0 + Qa + 18446))/1000 + (4*x15*x18*x21)/((x15 + 10)*(x21 + 1/5)) + (16*x15*x18*x22)/(25*(x15 + 10)*(x22 + 1/2)*(x21 + 1/5))),
            1/96*((x6*(Q0 + Qa + 18446))/1000 - x19/20 - (x19*(Q0 + Qa + 18446))/1000 + (x19*x21*x23)/(2*(x23 + 1)*(x21 + 2/5))),
            1/96*((3*x18)/125 + x19/250 + (x7*(Q0 + Qa + 18446))/1000 - (x20*(Q0 + Qa + 18446))/1000),
            1/96*((x8*(Q0 + Qa + 18446))/1000 - (x21*(Q0 + Qa + 18446))/1000 - (132*x15*x18*x21)/(67*(x15 + 10)*(x21 + 1/5)) - (433*x19*x21*x23)/(48*(x23 + 1)*(x21 + 2/5))),
            1/96*((x9*(Q0 + Qa + 18446))/1000 - (x22*(Q0 + Qa + 18446))/1000 + (25*x19*x21*x23)/(12*(x23 + 1)*(x21 + 2/5)) - (96*x15*x18*x22)/(871*(x15 + 10)*(x22 + 1/2)*(x21 + 1/5))),
            1/96*((x18*x24)/20 + (x10*(Q0 + Qa + 18446))/1000 - (x23*(Q0 + Qa + 18446))/1000 - (8*x15*x18*x21)/(25*(x15 + 10)*(x21 + 1/5)) - (637*x19*x21*x23)/(300*(x23 + 1)*(x21 + 2/5)) - (32*x15*x18*x22)/(625*(x15 + 10)*(x22 + 1/2)*(x21 + 1/5))),
            1/96*((x11*(Q0 + Qa + 18446))/1000 - (x18*x24)/20 - (x24*(Q0 + Qa + 18446))/1000 + (3*x25*(x21/(x21 + 1/5) + (4*x22)/(25*(x22 + 1/2)*(x21 + 1/5))))/(x17/x18 + 1/10)),
            1/96*((141*x18)/6250 + (47*x19)/12500 + (x12*(Q0 + Qa + 18446))/1000 - (x25*(Q0 + Qa + 18446))/1000 - (3*x25*(x21/(x21 + 1/5) + (4*x22)/(25*(x22 + 1/2)*(x21 + 1/5))))/(x17/x18 + 1/10)),
            1/96*((x18*x24)/280 + (x13*(Q0 + Qa + 18446))/1000 - (x26*(Q0 + Qa + 18446))/1000 - (4*x15*x18*x21)/(175*(x15 + 10)*(x21 + 1/5)) - (631*x19*x21*x23)/(2100*(x23 + 1)*(x21 + 2/5)) + (1518823277841921*x15*x18*x22)/(360287970189639680*(x15 + 10)*(x22 + 1/2)*(x21 + 1/5))),

            1/96*((x14*(Q0 + Qa + 18446))/1333 - (x27*(Q0 + Qa + 18446))/1333),
            1/96*((x15*(Q0 + Qa + 18446))/1333 - (x28*(Q0 + Qa + 18446))/1333 + (3*x30*(x34/(x34 + 1/5) + (4*x35)/(25*(x35 + 1/2)*(x34 + 1/5))))/(x30/x31 + 1/10) - (400*x28*x31*x34)/(67*(x28 + 10)*(x34 + 1/5)) - (64*x28*x31*x35)/(67*(x28 + 10)*(x35 + 1/2)*(x34 + 1/5))),
            1/96*((x16*(Q0 + Qa + 18446))/1333 - (x29*(Q0 + Qa + 18446))/1333),
            1/96*((69*x31)/250 + (23*x32)/500 + (x17*(Q0 + Qa + 18446))/1333 - (x30*(Q0 + Qa + 18446))/1333 - (3*x30*(x34/(x34 + 1/5) + (4*x35)/(25*(x35 + 1/2)*(x34 + 1/5))))/(x30/x31 + 1/10)),
            1/96*((x18*(Q0 + Qa + 18446))/1333 - (3*x31)/10 - (x31*(Q0 + Qa + 18446))/1333 + (4*x28*x31*x34)/((x28 + 10)*(x34 + 1/5)) + (16*x28*x31*x35)/(25*(x28 + 10)*(x35 + 1/2)*(x34 + 1/5))),
            1/96*((x19*(Q0 + Qa + 18446))/1333 - x32/20 - (x32*(Q0 + Qa + 18446))/1333 + (x32*x34*x36)/(2*(x36 + 1)*(x34 + 2/5))),
            1/96*((3*x31)/125 + x32/250 + (x20*(Q0 + Qa + 18446))/1333 - (x33*(Q0 + Qa + 18446))/1333),
            1/96*((x21*(Q0 + Qa + 18446))/1333 - 240*x34 - (x34*(Q0 + Qa + 18446))/1333 - (132*x28*x31*x34)/(67*(x28 + 10)*(x34 + 1/5)) - (433*x32*x34*x36)/(48*(x36 + 1)*(x34 + 2/5)) + 1920),
            1/96*((x22*(Q0 + Qa + 18446))/1333 - (x35*(Q0 + Qa + 18446))/1333 + (25*x32*x34*x36)/(12*(x36 + 1)*(x34 + 2/5)) - (96*x28*x31*x35)/(871*(x28 + 10)*(x35 + 1/2)*(x34 + 1/5))),
            1/96*((x31*x37)/20 + (x23*(Q0 + Qa + 18446))/1333 - (x36*(Q0 + Qa + 18446))/1333 - (8*x28*x31*x34)/(25*(x28 + 10)*(x34 + 1/5)) - (637*x32*x34*x36)/(300*(x36 + 1)*(x34 + 2/5)) - (32*x28*x31*x35)/(625*(x28 + 10)*(x35 + 1/2)*(x34 + 1/5))),
            1/96*((x24*(Q0 + Qa + 18446))/1333 - (x31*x37)/20 - (x37*(Q0 + Qa + 18446))/1333 + (3*x38*(x34/(x34 + 1/5) + (4*x35)/(25*(x35 + 1/2)*(x34 + 1/5))))/(x30/x31 + 1/10)),
            1/96*((141*x31)/6250 + (47*x32)/12500 + (x25*(Q0 + Qa + 18446))/1333 - (x38*(Q0 + Qa + 18446))/1333 - (3*x38*(x34/(x34 + 1/5) + (4*x35)/(25*(x35 + 1/2)*(x34 + 1/5))))/(x30/x31 + 1/10)),
            1/96*((x31*x37)/280 + (x26*(Q0 + Qa + 18446))/1333 - (x39*(Q0 + Qa + 18446))/1333 - (4*x28*x31*x34)/(175*(x28 + 10)*(x34 + 1/5)) - (631*x32*x34*x36)/(2100*(x36 + 1)*(x34 + 2/5)) + (1518823277841921*x28*x31*x35)/(360287970189639680*(x28 + 10)*(x35 + 1/2)*(x34 + 1/5))),

            1/96*((x27*(Q0 + Qa + 18446))/1333 - (x40*(Q0 + Qa + 18446))/1333),
            1/96*((x28*(Q0 + Qa + 18446))/1333 - (x41*(Q0 + Qa + 18446))/1333 + (3*x43*(x47/(x47 + 1/5) + (4*x48)/(25*(x48 + 1/2)*(x47 + 1/5))))/(x43/x44 + 1/10) - (400*x41*x44*x47)/(67*(x41 + 10)*(x47 + 1/5)) - (64*x41*x44*x48)/(67*(x41 + 10)*(x48 + 1/2)*(x47 + 1/5))),
            1/96*((x29*(Q0 + Qa + 18446))/1333 - (x42*(Q0 + Qa + 18446))/1333),
            1/96*((69*x44)/250 + (23*x45)/500 + (x30*(Q0 + Qa + 18446))/1333 - (x43*(Q0 + Qa + 18446))/1333 - (3*x43*(x47/(x47 + 1/5) + (4*x48)/(25*(x48 + 1/2)*(x47 + 1/5))))/(x43/x44 + 1/10)),
            1/96*((x31*(Q0 + Qa + 18446))/1333 - (3*x44)/10 - (x44*(Q0 + Qa + 18446))/1333 + (4*x41*x44*x47)/((x41 + 10)*(x47 + 1/5)) + (16*x41*x44*x48)/(25*(x41 + 10)*(x48 + 1/2)*(x47 + 1/5))),
            1/96*((x32*(Q0 + Qa + 18446))/1333 - x45/20 - (x45*(Q0 + Qa + 18446))/1333 + (x45*x47*x49)/(2*(x49 + 1)*(x47 + 2/5))),
            1/96*((3*x44)/125 + x45/250 + (x33*(Q0 + Qa + 18446))/1333 - (x46*(Q0 + Qa + 18446))/1333),
            1/96*((x34*(Q0 + Qa + 18446))/1333 - 240*x47 - (x47*(Q0 + Qa + 18446))/1333 - (132*x41*x44*x47)/(67*(x41 + 10)*(x47 + 1/5)) - (433*x45*x47*x49)/(48*(x49 + 1)*(x47 + 2/5)) + 1920),
            1/96*((x35*(Q0 + Qa + 18446))/1333 - (x48*(Q0 + Qa + 18446))/1333 + (25*x45*x47*x49)/(12*(x49 + 1)*(x47 + 2/5)) - (96*x41*x44*x48)/(871*(x41 + 10)*(x48 + 1/2)*(x47 + 1/5))),
            1/96*((x44*x50)/20 + (x36*(Q0 + Qa + 18446))/1333 - (x49*(Q0 + Qa + 18446))/1333 - (8*x41*x44*x47)/(25*(x41 + 10)*(x47 + 1/5)) - (637*x45*x47*x49)/(300*(x49 + 1)*(x47 + 2/5)) - (32*x41*x44*x48)/(625*(x41 + 10)*(x48 + 1/2)*(x47 + 1/5))),
            1/96*((x37*(Q0 + Qa + 18446))/1333 - (x44*x50)/20 - (x50*(Q0 + Qa + 18446))/1333 + (3*x51*(x47/(x47 + 1/5) + (4*x48)/(25*(x48 + 1/2)*(x47 + 1/5))))/(x43/x44 + 1/10)),
            1/96*((141*x44)/6250 + (47*x45)/12500 + (x38*(Q0 + Qa + 18446))/1333 - (x51*(Q0 + Qa + 18446))/1333 - (3*x51*(x47/(x47 + 1/5) + (4*x48)/(25*(x48 + 1/2)*(x47 + 1/5))))/(x43/x44 + 1/10)),
            1/96*((x44*x50)/280 + (x39*(Q0 + Qa + 18446))/1333 - (x52*(Q0 + Qa + 18446))/1333 - (4*x41*x44*x47)/(175*(x41 + 10)*(x47 + 1/5)) - (631*x45*x47*x49)/(2100*(x49 + 1)*(x47 + 2/5)) + (1518823277841921*x41*x44*x48)/(360287970189639680*(x41 + 10)*(x48 + 1/2)*(x47 + 1/5))),

            1/96*((x40*(Q0 + Qa + 18446))/1333 - (x53*(Q0 + Qa + 18446))/1333),
            1/96*((x41*(Q0 + Qa + 18446))/1333 - (x54*(Q0 + Qa + 18446))/1333 + (3*x56*(x60/(x60 + 1/5) + (4*x61)/(25*(x61 + 1/2)*(x60 + 1/5))))/(x56/x57 + 1/10) - (400*x54*x57*x60)/(67*(x54 + 10)*(x60 + 1/5)) - (64*x54*x57*x61)/(67*(x54 + 10)*(x61 + 1/2)*(x60 + 1/5))),
            1/96*((x42*(Q0 + Qa + 18446))/1333 - (x55*(Q0 + Qa + 18446))/1333),
            1/96*((69*x57)/250 + (23*x58)/500 + (x43*(Q0 + Qa + 18446))/1333 - (x56*(Q0 + Qa + 18446))/1333 - (3*x56*(x60/(x60 + 1/5) + (4*x61)/(25*(x61 + 1/2)*(x60 + 1/5))))/(x56/x57 + 1/10)),
            1/96*((x44*(Q0 + Qa + 18446))/1333 - (3*x57)/10 - (x57*(Q0 + Qa + 18446))/1333 + (4*x54*x57*x60)/((x54 + 10)*(x60 + 1/5)) + (16*x54*x57*x61)/(25*(x54 + 10)*(x61 + 1/2)*(x60 + 1/5))),
            1/96*((x45*(Q0 + Qa + 18446))/1333 - x58/20 - (x58*(Q0 + Qa + 18446))/1333 + (x58*x60*x62)/(2*(x62 + 1)*(x60 + 2/5))),
            1/96*((3*x57)/125 + x58/250 + (x46*(Q0 + Qa + 18446))/1333 - (x59*(Q0 + Qa + 18446))/1333),
            1/96*((x47*(Q0 + Qa + 18446))/1333 - KLa5*(x60 - 8) - (x60*(Q0 + Qa + 18446))/1333 - (132*x54*x57*x60)/(67*(x54 + 10)*(x60 + 1/5)) - (433*x58*x60*x62)/(48*(x62 + 1)*(x60 + 2/5))),
            1/96*((x48*(Q0 + Qa + 18446))/1333 - (x61*(Q0 + Qa + 18446))/1333 + (25*x58*x60*x62)/(12*(x62 + 1)*(x60 + 2/5)) - (96*x54*x57*x61)/(871*(x54 + 10)*(x61 + 1/2)*(x60 + 1/5))),
            1/96*((x57*x63)/20 + (x49*(Q0 + Qa + 18446))/1333 - (x62*(Q0 + Qa + 18446))/1333 - (8*x54*x57*x60)/(25*(x54 + 10)*(x60 + 1/5)) - (637*x58*x60*x62)/(300*(x62 + 1)*(x60 + 2/5)) - (32*x54*x57*x61)/(625*(x54 + 10)*(x61 + 1/2)*(x60 + 1/5))),
            1/96*((x50*(Q0 + Qa + 18446))/1333 - (x57*x63)/20 - (x63*(Q0 + Qa + 18446))/1333 + (3*x64*(x60/(x60 + 1/5) + (4*x61)/(25*(x61 + 1/2)*(x60 + 1/5))))/(x56/x57 + 1/10)),
            1/96*((141*x57)/6250 + (47*x58)/12500 + (x51*(Q0 + Qa + 18446))/1333 - (x64*(Q0 + Qa + 18446))/1333 - (3*x64*(x60/(x60 + 1/5) + (4*x61)/(25*(x61 + 1/2)*(x60 + 1/5))))/(x56/x57 + 1/10)),
            1/96*((x57*x63)/280 + (x52*(Q0 + Qa + 18446))/1333 - (x65*(Q0 + Qa + 18446))/1333 - (4*x54*x57*x60)/(175*(x54 + 10)*(x60 + 1/5)) - (631*x58*x60*x62)/(2100*(x62 + 1)*(x60 + 2/5)) + (1518823277841921*x54*x57*x61)/(360287970189639680*(x54 + 10)*(x61 + 1/2)*(x60 + 1/5))),
 
#            1/96*((6277*x67)/200 - (6277*x66)/200 + (5*x67*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x67)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x67)/50000)))/2),
#            1/96*((6277*x68)/200 - (6277*x67)/200 - (5*x67*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x67)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x67)/50000)))/2 + (5*x68*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x68)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x68)/50000)))/2),
#            1/96*((6277*x69)/200 - (6277*x68)/200),
#            1/96*((6277*x70)/200 - (6277*x69)/200 - (5*x68*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x68)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x68)/50000)))/2 + (5*x70*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x70)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x70)/50000)))/2),
# SI SS XI XS XBH XBA XP SO SNO SNH SND XND SALK
            1/96*((x53*(Q0 + 18446))/600 - (5*x72*(Q0/1500 + 9223/750))/2), 
            1/96*((x54*(Q0 + 18446))/600 - (5*x73*(Q0/1500 + 9223/750))/2),
            
            1/96*((6277*x67)/200 - (6277*x66)/200),
            1/96*((5*x68*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x68)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x68)/50000)))/2 - (5*x66*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x66)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x66)/50000)))/2 + ((Q0 + 18446)*((3*x55)/4 + (3*x56)/4 + (3*x57)/4 + (3*x58)/4 + (3*x59)/4))/600 - (5*x67*(Q0/1500 + 9223/750))/2), 
            1/96*((5*(Q0/1500 - 77/300)*(x67 - x68))/2 - (5*x68*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x68)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x68)/50000)))/2 + (5*x69*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x69)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x69)/50000)))/2),
            1/96*((5*(Q0/1500 - 77/300)*(x68 - x69))/2 - (5*x69*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x69)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x69)/50000)))/2 + (5*x70*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x70)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x70)/50000)))/2),
            1/96*((5*(Q0/1500 - 77/300)*(x69 - x70))/2 - (5*x70*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x70)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x70)/50000)))/2 + (5*x71*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x71)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x71)/50000)))/2),

            1/96*((x60*(Q0 + 18446))/600 - (5*x74*(Q0/1500 + 9223/750))/2),
            1/96*((x61*(Q0 + 18446))/600 - (5*x75*(Q0/1500 + 9223/750))/2),
            1/96*((x62*(Q0 + 18446))/600 - (5*x76*(Q0/1500 + 9223/750))/2),
            1/96*((x63*(Q0 + 18446))/600 - (5*x77*(Q0/1500 + 9223/750))/2),
            
            1/96*((5*(Q0/1500 - 77/300)*(x70 - x71))/2 - (5*x71*(474* np.exp((1539*x55)/1562500000 + (1539*x56)/1562500000 + (1539*x57)/1562500000 + (1539*x58)/1562500000 + (1539*x59)/1562500000 - (9*x71)/15625) - 474* np.exp((24453*x55)/5000000000 + (24453*x56)/5000000000 + (24453*x57)/5000000000 + (24453*x58)/5000000000 + (24453*x59)/5000000000 - (143*x71)/50000)))/2),
            
            1/96*((x65*(Q0 + 18446))/600 - (5*x78*(Q0/1500 + 9223/750))/2),]
    
    return np.array(dx_scaledt)
#Ny = 13
#def measurement_wwtp(x):
#    C_extend = np.zeros([Ny,Nx])
##    C_extend[0,20] = np.array([1])     #SO2,3,4,5
##    C_extend[1,33] = np.array([1])
#    C_extend[0,46] = np.array([1])
#    C_extend[1,59] = np.array([1])
#    C_extend[2,[13,14,15,16,17,18]] = np.array([1,1,1,1,1,1])  # COD2, COD3
#    C_extend[3,[26,27,28,29,30,31]] = np.array([1,1,1,1,1,1])
#    C_extend[4,[65,66,67,68,71,72]] = np.array([1,1,1,1,1,1])  # COD6
#    C_extend[5,[13,14]] = np.array([1,1]) # CDOf2, CODf3
#    C_extend[6,[26,27]] = np.array([1,1])
##    C_extend[9,25] = np.array([1])        # SALK2, SALK3, SALK6
#    C_extend[7,38] = np.array([1])     
#    C_extend[8,77] = np.array([1])
#    C_extend[9,[15,16,17,18,19,24]] = np.array([1,1,1,1,1,1]) # TSS2
#    C_extend[10,[28,29,30,31,32,37]] = np.array([1,1,1,1,1,1]) # TSS3
#    C_extend[11,[74]] = np.array([1])      
#    C_extend[12,[75]] = np.array([1])
#    y = np.dot(C_extend,x)
#    return np.array(y)

